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We develop an effective theory of pulse propagation in a nonlinear and 
disordered medium. The theory is formulated in terms of a nonlinear dif- 
fusion equation. Despite its apparent simplicity this equation describes 
novel phenomena which we refer to as "locked explosion" and "diffu- 
sive" collapse. The equation can be applied to such distinct physical 
systems as laser beams propagating in disordered photonic crystals or 
Bose-Einstein condensates expanding in a disordered environment. 

1.1. Introduction 

In recent years, novel experimental techniques made possible first obser- 
vations of wave-packets evolving in the presence of random scatterers and 
nonlinearities. In a number of optical experiments, a laser beam was sent 
into a nonlinear optical medium with a random refractive index, and the 
beam profile in the transverse direction(s) was monitored on the opposite 
side of the sample, 1 2 for an illustration see Fig. |l.l[ a). In a second class 
of experiments, atoms forming a Bose-Einstein condensate were released 
from a trap and subjected to a disorder potential during the expansionpH^ 
see Fig. |l.l[ b). The experiments were inspired by the idea that in these 
setups, unlike for transport experiments in electronic systems, one can vi- 
sualize the phenomenon of Anderson localization, whereby a wave-packet 
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(a) (b) 

Fig. 1.1. (a) In the experiment of Ref. [Ila laser beam (in red) is sent onto a disordered 
nonlinear photonic crystal. Inside the crystal, the evolution of the smooth envelope of 
the electric field in the transverse (xy) plane is approximately described by the nonlinear 
Schrodinger equation (NLSE). In this equation, the z-coordinate plays the role of "time". 
The intensity profile of the outgoing beam is measured at the opposite side of the crystal. 
Adapted from Schwartz et al., Ref. [T] (b) In the experiment of Ref. [7] an atomic Bose- 
Einstein condensate is released from a small trap (displayed in white) and subsequently 
expands along a one-dimensional channel, formed by an external confining potential (in 
red). During the expansion, the atoms are subjected to a disorder potential (in blue). 
After a while, the motion comes to a halt and the condensate becomes localized. The 
experiment was performed at small densities, so that the nonlinear term in the Gross- 
Pitaevskii equation can be considered as small. Adapted from Billy et al., Ref. [7] 



or quantum particle is confined within a finite volume as a result of mul- 
tiple scattering on a random potential (Figs. 1.2|1.3|1.4 ). The evolution 
of the injected wave-packet in both experiments can be described by the 
non-linear Schrodinger equation (NLSE), in the context of atomic Bose- 
Einstein condensates referred to as the Gross-Pitaevskii Equation (GPE). 
This equation differs from the linear Schrodinger equation by an additional 
cubic term and is used as a paradigmatic description for nonlinear waves. 
The nonlinearity is a consequence of interactions between particles in the 
case of atomic condensates and of a change in the refractive index in re- 
sponse to the electric field (Kerr effect) in the case of laser beams. 

Motivated by these experiments we derive, starting from the 
GPE/NLSE, a kinetic equation that describes the evolution of an injected 
wave-packet in a weakly disordered nonlinear medium in two dimensions. 
Analysis of this equation reveals a rather nontrivial picture: Irrespective 
of the sign of the nonlinearity the mean square radius of the wave-packet 
changes linearly in time, dt (r 2 ) oc E tot , where E tot is the total energy 
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of the wave-packet. For a repulsive nonlinearity the initial change of the 
profile displays features of an explosion, although the overall size of the 
wave-packet is growing slowly as in ordinary diffusion. For an attractive 
nonlinearity, the radius can either grow or decrease, depending on the sign 
of E to t- I n particular, for E tot < we predict a slow "diffusive" collapse as 
the radius of the wave-packet shrinks towards zero. 

In this paper we will mostly use the language related to the GPE, but 
also indicate below how to translate to a language more suitable for optical 
experiments. The GPE describes the evolution of a wave- function \I/: 9 

ifc¥(r,t) = --^V 2 *(r,t) + u(r)*(r,t) + A|*(r,t)| 2 *(r,t), (1.1) 
2m 

where we set ft = 1. For positive (negative) A this equation contains a 
repulsive (attractive) self-consistent potential A|\I/(r, £)| 2 . This corresponds 
to a nonlinearity of the de-focusing (self-focusing) type. The static disorder 
potential u(r) is the source of randomness in the equation. For simplicity we 
choose for our calculation a Gaussian white noise potential with correlation 
function (u(r)u(r')} = S(r — r / )/(mr) [for a discussion of averaging for 
speckle potentials see, e.g., Ref. 10 . The angular brackets denote averaging 
over disorder configurations and r is the scattering time. 

The NLSE used in optics is derived in the so-called paraxial approxima- 




Fig. 1.2. On the left hand side, the profile of the condensate at the final stage of the 
experiment of Ref. [7] is displayed on a logarithmic scale. The tails of the wave function 
decay exponentially. This is interpreted as evidence for Anderson localization. On the 
right hand side the rms width of the condensate is plotted as a function of time. In the 
presence of a disorder potential (red line), the width approaches a constant value, while 
in the absence of disorder (green line), the condensate expands ballistically. Adapted 
from Billy et al., Ref. [7] 
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Fig. 1.3. The intensity distribution of the outgoing laser beam in the experiment of 
Ref. ^ is shown in a regime for which the nonlinearity is negligible. The disorder level 
increases from left to right, starting from the clean lattice. All distributions are aver- 
aged over different disorder configurations. The white lines display the logarithm of the 
intensity as a function of the transverse coordinate. The authors of Ref. 1 interpret the 
distribution in the middle as the result of diffusion, and the distribution on the right as 
a signature of Anderson localization. Courtesy of Schwartz et al., Ref. 1 

tionpl and describes the evolution of the smooth envelope of the electric 
field. The main propagation direction of the laser beam, say the z-direction, 
plays the role of time in the NLSE, see Fig.[DJa). In this sense, the disor- 
der potential which results from random variations of the refractive index 
is static when it is ^-independent. For example, the two-dimensional (2d) 
transverse evolution of a pulse is studied in a 3d sample. The intensity of 
the beam is proportional to |\I/(r, z)\ 2 . In the NLSE, the mass m in the 
GPE is replaced by the wave vector k = cj/c, where uo is the frequency of 
the carrier wave and c the velocity of light in the medium. 

For a condensate released from a confining harmonic oscillator poten- 
tial, as is typical for experiments on cold atomic gases, the GPE without 
disorder can be solved exact lyllHlHI During an initial stage the potential 
energy originating from the nonlinearity is almost entirely converted into 
kinetic energy. This period of violent acceleration is followed by a sec- 
ond stage, during which the nonlinearity is no longer essential. Expansion 
in the presence of disorder in the two-dimensional case was recently ad- 
dressed in reference.^ In this paper it has been assumed that for repulsive 
nonlinearity an initial ballistic stage is not affected by disorder, while the 
subsequent diffusive expansion is not affected by the nonlinearity, thereby 
separating the two effects. In contrast, we are interested here in the in- 
terplay of disorder and nonlinearity, both attractive and repulsive. This 
is especially interesting in 2d, as it is known that for linear wave propa- 
gation and weak disorder there is an extended diffusive regime preceding 
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Fig. 1.4. The average width of the outgoing laser beam as a function of the disorder 
level in the experiment of Ref. [T] is compared for the linear case and in the presence of a 
self-focusing nonlinearity (left panel). The right panel compares (the logarithm of) the 
averaged intensity profile for a fixed disorder strength as a function of the transverse 
coordinate. The parameter a is a dimensionless measure for the strength of the non- 
linearity. The main result is that the self-focusing nonlinearity promotes localization. 
Courtesy of Schwartz et al., Ref. ^| 

localization on (exponentially) large length scales. It is this regime that we 
address 

Since the system is far out of equilibrium, we choose to work with a 
kinetic equation. The derivation of the kinetic equation proceeds as follows. 
We use methods of classical statistical field theory to derive a functional 
integral expression for the disorder averaged density U£LLZI The formalism 
involves a doubling of the degrees of freedom, similar to the Keldysh or 
closed-time-path approaches for quantum systems where two fields are 
introduced on forward and backward time-contours. Instead of averaging 
over a statistical ensemble in the initial state, we assume that the wave- 
function at the initial time is known. Averaging is performed over disorder 
configurations. Scattering on impurities is included on the level of the 
self-consistent Born approximation. While interference (weak localization) 
corrections are not covered by this approximation, it allows for a consistent 
description of diffusion in the presence of nonlinearity. The nonlinearity 
is treated by introducing a self-consistent potential #(r, t). In this way it 
is possible to include interaction effects in a non-perturbative way, which 
is crucial for the problem at hand. To obtain the kinetic equation for 
the density in the diffusive limit, we assume that the initial wave- function 
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sets a momentum scale po characterizing the main part of the momentum 
distribution, so that the weak disorder condition pol ^> 1 is fulfilled, where 
I = Por/m is the mean free path. We further assume that the density 
varies smoothly on scales of Z, in particular that the size of the condensate 
is much larger than the mean free path. Both of these conditions can be 
met simultaneously. The phase of which is related to the momentum, 
may change rapidly, while the amplitude, which determines the density, 
may vary smoothly. Even if the density does not satisfy the smoothness 
condition initially, it is natural to expect that in the case of an expansion 
it will become sufficiently smooth after some timeP^ The derivation of the 
kinetic equation will be presented elsewhere.^ 

Starting from Eq. the outlined steps lead to the following kinetic 

equation in the diffusive regime 

dtn(r, t, e) - V(L> £ _tf Vn(r, t, e)) + d t #(r, t)d e n(r, t, e) 

= 5(t)F(e-#(r,0),r) (1.2) 



where 



F(s, r) = J ±0-£ F(p, q) exp(iqr) 2*5{e - e p ), (1.3) 

and F(p,q) = ^o(p + q/2)^o(P — is determined by the initial wave 
function = p 2 /(2m) is the kinetic energy, and D £ = er/m the 

diffusion coefficient. The equation should be supplemented with the self- 
consistency relation for the potential $(r,£) = 2An(r,£), where n(r, t) = 
J ds/(27r) n(r,t,e). Despite its apparent simplicity it is a rather compli- 
cated nonlinear integro-differential equation. The equation effectively sums 
an infinite series of diagrams of the type shown in Fig. |1.5| It is a peculiarity 
of the perturbation theory for a classical field equation such as the GPE that 
no closed loops arise™ making it quite distinct from the related problem in 
interacting electron systems. The relation between certain blocks appear- 
ing in diagrammatic perturbation theory and the corresponding terms in 



the kinetic equation is visualized in Fig. 1.6 

The physics described by this equation is essentially classical. Imagine 
first that the potential $ does not depend on time. Consider now a particle 
diffusing with total energy e on the background of a smoothly varying 



potential see Fig. 1.7 for illustration. If scattering events are frequent 
enough, the diffusion coefficient is determined by the kinetic energy e p = 
e — d that varies locally in space. If the potential additionally varies in time, 
the particle may change its total energy. If on the other hand the potential 
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Fig. 1.5. Graphical illustration of the kinetic equation, Eq. ( |1.2| ). The injection pro- 
cess takes place on the left. Solid lines are disorder averaged particle propagators. Lad- 
ders graphically represent density diffusion. The particular way of disorder averaging 
is justified for EkinT ^> 1; ^kin is defined below Eq. (1.6) . The wavy lines account for 
the nonlinearity in Eq. ( |l.l| ). The block magnified in the inset gives rise to the terms 
V(D^Vn(r,t,e)) and dt$(r,t)d £ n(r,t, e) in the kinetic equation, Eq. 2. 



d t h{r, t, e) - V( J D £ _ 1? Vn(r, t, e)) + d t $(r, t)d £ h(r, t, e) = S(t) F(e - <&(r, 0), r) 



I \ 




Fig. 1.6. This figure illustrates the correspondence between different terms in the kinetic 
equation and a certain block in perturbation theory. 

depends on time only, the kinetic energy does not change. Therefore, it 
is expected that a purely time-dependent potential has no effect on the 
density. This observation is related to the fact that in the original GPE a 
purely time dependent potential V(t) may be removed by a suitably chosen 
gauge-transformation, 

tf(r,t) -> tf(r,t) exp ^ dt'V(t')^j , (1.4) 
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that does not affect the density |\I/(r, t)\ 2 . Indeed, we can make this 
point obvious in Eq. (1.2) by shifting the energy variable so that it will 
correspond to the kinetic energy instead of the total energy, n(r, e, t) = 
n(r, s + $(r, £),£). Expressed in the new coordinates the equation reads 

dtn(r, e, t) - [V - W r , A] D e [V - W P , t e ] n(r, e, t) = S(t) F(e, r) (1.5) 




(a) 



Fig. 1.7. The density distribution n(r, t) creates its own self-consistent potential #(r, t). 
We illustrate the attractive (focusing) case. In a classical analogy, the variable e in Eq. |1.2| 
can be interpreted as the total energy and •& as the potential energy of a diffusing particle. 
Correspondingly, the diffusion coefficient is determined by the kinetic energy e p = e — ti. 

An equation for the density n(r,t) can be obtained by integrating 



Eq. ([L5j) in e 

d t n(r, t) - —V 2 (e(r, t) + An 2 (r, t)) = n(r, 0), (1.6) 



m 

where 



f de 

e(r,t)= / — en(r,t,e) = e fein (r,t)n(r,t). (1.7) 

It can be written in the compact form d t n — \7 2 (D e ffn) = S(t)n when 
defining an effective space and time-dependent diffusion coefficient D e ff = 
(£kin + Xn)r/m. The apparent simplicity of these equations is however 
deceiving. They are not closed equations for the density evolution, since 
the kinetic energy e depends on the nonlinearity and needs to be deter- 



mined separately via Eq. (1.5). Nevertheless, we arrive at the conceptually 
important result that in the diffusive regime the nonlinearity effectively in- 
troduces a density dependence of the diffusion coefficient. It seems clear 
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Fig. 1.8. This figure compares two different linear diffusion processes, one with a fixed 
diffusion coefficient (gray) and one with energy-dependent diffusion coefficient D £ and an 
initial energy distribution of the form Cexp(— s/sq), where C and £q are constants. On 
the left hand side the initial density distribution is shown, which is chosen to be identical 
for both processes. After a while, the density that corresponds to energy dependent 
diffusion is larger both for small distances and for large distances from the origin when 
compared to the diffusion process with a fixed diffusion constant D £Q . This reflects the 
fact that both energies e <C £o and e ^> £q are present in the initial energy distribution, 
leading to slow and fast diffusion respectively. 



that a closed form solution of the nonlinear equations for arbitrary initial 
conditions cannot be found. In order to make progress we will rely on two 
approaches: the use of conservation laws and the study of solvable limiting 
cases. When combined, they will enable us to arrive at a qualitative picture 
both for repulsive and attractive nonlinearity. 

First we briefly discuss the linear case, $ = 0. In the absence of non- 
linearity, n(r,£,t) evolves independently for each energy e. In this limit, 
Eq. (1.5) has the obvious solution 

n(r, e, t) = ^ J dr, e~^/^F{e, n). (1.8) 

Here, for each energy e diffusion is determined by the corresponding dif- 
fusion coefficient D £1 and should be weighted according to the energy dis- 
tribution in the injected wave-packet. The linear case was discussed in 



Ref. 14 Starting from a broad energy distribution centered around some 
So, both tails and central part of the density in the long time limit are 
more pronounced compared to diffusion at fixed energy eq, because the 
tails are determined by energies e > while the central part is dominated 



by energies e < see Fig. 1. 



Next we turn to the nonlinear case. We will make use of the conservation 
laws for particle number and energy. By integrating Eq. ( |1.6| ) over r, we 
obtain that the particle number (or normalization) J dr n(r, t) = N is fixed 
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in time. An equation for e(r, t) can be derived by first multiplying Eq. ( 1.5 ) 
by e before integrating in this variable. Then by combining the equation 
for n(r, t) with the equation for e(r, t) we find that the energy 



J dr (e(r,t) + \n 2 (r,t)) (1.9) 



is constant in time. Remarkably, this conservation law completely de- 
termines the time evolution of the mean radius squared of the wave- 
packet, (r 2 ) = J dr r 2 n(r,t)/N. Indeed, multiplying Eq. (1.6) by r 2 



and subsequently integrating in r one obtains that dt (r 2 ) = 4D £tot , where 
Stot — Etot/N. The linear dependence of the mean square radius on time 
during the whole evolution is guarded by energy conservation. This is one 
of the central results of this paper. When compared to the linear case, the 
effective diffusion coefficient D £tot is reduced for attractive and enhanced 
for repulsive nonlinearities. 

In the following we discuss more specifically the repulsive and attractive 
cases. For the repulsive nonlinear case it is instructive to consider a situ- 



ation in which the second term on the RHS of Eq. (1.6) dominates. The 
equation d t n = V 2 n 2 , which one obtains after simple rescaling, is an exam- 
ple of the famous porous medium equation (PME) PH For the 2d case the 
solution describing the evolution of a delta- function pulse MS(r) is given 
by n(r,£) = (C - r 2 /(Wi 1 / 2 ))/^ 2 , where C 2 = M/^) 22 ^ This solution 
is often referred to as Barenblatt's solution. It conserves the normalization 
J dr n(r, t) = M but, unlike ordinary diffusion, it is nonzero only in a finite 
region of space, see Fig. |1.9| The special importance of Barenblatt's solu- 
tion in the theory of the PME is related to the fact that, roughly speaking, 
any solution starting from a sufficiently benign initial pulse with weight 
M is eventually well- approximated by Barenblatt's solution with the same 
weight 

For Barenblatt's solution, the mean radius squared evolves as (r 2 ) cx 
t x l 2 and the density at r = drops as n(0,t) cx t~ x l 2 . At short times this 
solution describes a much faster "explosive" evolution than the source- type 
solution of the diffusion equation, for which (r 2 ) cx t and 77,(0, t) cx t -1 , 
see Fig. |1.1Q| At first sight there seems to be a contradiction. If one 
injects a bell-shaped pulse with a large potential energy, it appears that the 
potential part of the effective diffusion coefficient D e ff = {euin + Xn)r/m 
dominates. Therefore, naively, one would assume that the initial evolution 
is "explosive", while our exact result (r 2 ) = Tq + AD £tot t rules out this 
possibility. This puzzle can be resolved in the following way. The explosion 
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Fig. 1.9. This figure shows the time-evolution of the Barenblatt-solution of the porous 
medium equation. It has the shape of an inverted parabola. It is worth noting that the 
wave-front has a finite spatial derivative at the boundary of the distribution. For this 
solution, the mean radius squared evolves rapidly as (r 2 ) oc t 1 / 2 and correspondingly 
the density at r = drops as n(0, t) oc t -1 / 2 . 
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Fig. 1.10. This figure compares a conventional diffusion process (black) to the time 
evolution described by the porous medium equation (gray) in two spatial dimensions. 
The same initial distribution is chosen for both processes and displayed on the left hand 
side. One can see, that the density at the center drops much faster for the solution of 
the porous medium equation. Simultaneously, the mean radius grows much faster, so 
that the total density is conserved. 



takes place only in the central part of the density distribution, which has 
only a small weight when calculating (r 2 ). Right in the center, for r = 0, 



Eq. (1.6) can be written as 

d t n = — \V 2 e + 2AnV 2 nl (1.10) 
rn 1 J 

for t > 0; we consider here a rotationally symmetric distribution with 
Vn(0,£) = and V 2 n(0,0) < 0. For sufficiently large An, the potential 
part is dominant and leads to a fast initial decrease of the density before 
either AnV 2 n becomes small or V 2 £ becomes positive as a consequence of 
the outward- flow of the kinetic energy. Away from the center, where the 



density and correspondingly the term 2AnV 2 n are small, Eq. (1.6) takes 
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the form 



d t n * — \V 2 e + 2A(Vn) 2 l (1.11) 
m 

for t > 0. For the PME it is the second term that determines the prop- 



agation of the boundary. For Eq. (1.6), however, the large kinetic energy 
outside the center leads to V 2 £ < for intermediate distances, and this 
prevents the term 2A(Vn) 2 from dominating. It is therefore an inversion of 
the distribution of kinetic energy compared to that of the density that does 
not allow for an explosive expansion and leads to a linear dependence of 
(r 2 ) on t. A sketch of a typical density evolution expected for this "locked 
explosion" is presented in the first line of Fig. 1.11| 



E tot >0 



^1 



Fig. 1.11. Time-evolution of the disorder averaged density n sketched qualitatively for 
the three relevant cases. According to the relation dt (r 2 ) = EtotT/m, the radius grows 
for Etot > and shrinks for Etot < 0. The condition Etot > can be realized for 
positive (first line) or negative potential $ = An (second line), while Etot < can be 
realized only for negative $ (third line). In the repulsive case a rapid drop of the density 
distribution is expected in the center for large # = An, while (r 2 ) grows only linearly 
in t ("locked explosion"). For Etot < a "diffusive" collapse is expected. In the cases 
with < a fragmentation of the cloud may occur. The evolution is determined by the 
initial density and energy distributions. 

We now discuss general features of wave-packet dynamics in the dis- 
ordered and nonlinear medium! 2 -^ For an expanding wave-packet, i.e. 
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Etot > 0, the overall potential energy related to the nonlinearity is con- 
verted into kinetic energy. As a result, the total kinetic energy increases 
in the repulsive case and decreases in the attractive case. Correspondingly, 
during the course of the expansion localization effects can be expected to be 
weakened for repulsive nonlinearity and enhanced for attractive nonlinear- 
ity. In particular, for an attractive (self-focusing) nonlinearity the slowing 
down and eventual localization of the injected pulse (not considered here) 
occurs at smaller distances than in the linear case as observed in the experi- 
ment PI [Regarding the role of interference effects for localization three main 
factors need to be accounted for. As long as the evolution does not come 
to a halt, the time-dependent potential leads to dephasing, which weakens 
localization effects. When the wave-packet becomes broader, longer paths 
become available for interference and at the same time the kinetic energy 
decreases, whereby sr decreases. The latter two effects support localiza- 
tion.] 

The attractive case is richer than the repulsive one (see Fig. 2), because 
the total energy may also be negative, E to t < 0. Then the mean radius 
squared would become equal to zero after a finite time. This corresponds 
to a celebrated phenomenon in nonlinear physics, the collapse EMU Here 
it is realized for the diffusive system. To the best to our knowledge, this 
"diffusive" collapse has not been discussed in the literature. Since our 
reasoning is based on a diffusive kinetic equation and thus assumes frequent 
scattering, the linear decrease of (r 2 ) only holds as long as the radius of 
the cloud exceeds the mean free path. [In the clean case the virial theorem 
for the NLSE in 2cl 25 l 26 l states that the second time derivative of (r 2 ) is 
proportional to the total energy. In this article we describe diffusive motion 
and correspondingly obtain a different time dependence for the size of the 
cloud.] 

Even for E to t > the collapse can play a role when the nonlinearity is 
attractive, if part of the cloud has a negative energy, while the remaining 
part expands. As a result one can expect a fragmentation of the cloud. 
If a part of the cloud with a positive but small energy lags behind, this 
fragment may have a strong tendency to localize. One may expect that 
this kind of localized or collapsing fragment generically remains from an 
expanding cloud with E tot > but attractive nonlinearity. 

To conclude, the nonlinear diffusion equation discussed in this paper 
contains rich physics that invites further theoretical and experimental in- 
vestigations. 
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